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Abstract 

A variational data assimilation technique was used to estimate 
optimal discretization of interpolation operators and derivatives in 
the nodes adjacent to the rigid boundary. Assimilation of artificially 
generated observational data in the shallow-water model in a square 
box and assimilation of real observations in the model of the Black 
sea are discussed. It is shown in both experiments that controlling 
the discretization of operators near a rigid boundary can bring the 
model solution closer to observations as in the assimilation window 
and beyond the window. This type of control allows also to improve 
climatic variability of the model. 

Keywords: Variational Data Assimilation; Boundary conditions; 
Shallow water model; Black sea model. 

1 Introduction 

Variational data assimilation technique, first proposed in [1], [2], is based on 
the optimal control methods [3] and perturbations theory This technique 
allows us to retrieve an optimal data for a given model from heterogeneous 
observation fields. Since the early 1990's, many mathematical and geophys- 
ical teams are involved in the development of the data assimilation strategy. 



*INRIA, projet MOISE, Laboratoire Jean Kuntzmann, BP 53, 38041 Grenoble Cedex 
9, France 



1 



One can cite many papers devoted to this problem, as in the domain of devel- 
opment of different methods for the data assimilation and also in the domain 
of its applications to the atmosphere and oceans. 

In the beginning, data assimilation methods were intended to identify and 
reconstruct an optimal initial state for the model. However, the idea that 
other model's parameters should also be identified by data assimilation has 
also been studied and discussed in numerous papers. One can cite several 
examples of using data assimilation to identify the bottom topography of 
simple models ([5], [6], [7]), to control open boundary conditions in coastal 
and regional models ([8], [9], [lO], [TJ) and to determine other parameters 
of a model ([12], [13], IH) • 

Together with these parameters, it seems to be interesting also to control 
boundary conditions on rigid boundaries. Although the boundary configura- 
tion of the ocean is steady and can be measured with much better accuracy 
than the model's initial state, it is not obvious how to represent it on the 
model's grid because of limited resolution. The coastal line of continents 
possesses a very fine structure and can only be roughly approximated by the 
model's grid. Consequently, boundary conditions are defined at the model 
grid's points which are different from the coast. Even the most evident im- 
permeability condition being placed at a wrong point may lead to some error 
in the model's solution. If it can improve the model solution, we may accept 
the flux can cross the boundary in places where the boundary is in water, 
prescribing some integral properties on the flux. 

Ocean models frequently include strong and thin boundary currents with 
intense velocity gradients. In this case, particular attention must be paid 
to the approximation of the boundary layer because a wrong representation 
of these currents may be responsible for drastic deformations of the global 
solution (see, for example [15] )■ This may lead us to control the discretization 
of the model's operators in the whole boundary layer rather than only at 
nodes directly adjacent to the boundary. 

Alternative method that is frequently used in geophysical models consists 
in the grid refinement in zones where boundary layers might occur. However, 
this implies additional computational cost on each model run. Variational 
data assimilation may help us to determine the parametrization of the bound- 
ary layer once for all model runs and to save the computer time. 

Boundary conditions on the rigid boundaries have been controlled by data 
assimilation for heat equation (see, for example, [16], [IZ]), but this control 
concerns the linear parabolic diffusion operator rather than hyperbolic trans- 
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port and advection operators that are more important in geophysical models. 

Studies on the possibility to control boundary conditions on rigid bound- 
aries for equations containing hyperbolic operators can be found in [18] on 
the example of non-linear balance equation, in [19] on the example of the 
wave equation, in [20] on the example of a linear shallow water model in a 
square box and in [21] and [22] on the example of the Burgers equation. The 
principal possibility to improve the model's solution controlling its boundary 
values are shown in all these papers. However, as it has been noted in [2T] . 
particular attention must be paid to the discretization process which must 
respect several rules. It is the discretization of the model's operators that 
takes into account the set of boundary conditions and introduces them into 
the model. Consequently, instead of controlling boundary conditions them- 
selves, there has been proposed in [20] to identify optimal discretization of 
differential operators at points adjacent to boundaries. This allows us to 
control directly the way the boundary conditions influence the model and to 
control boundary parameters in a more general way. Boundary conditions 
participate in discretized operators, but considering the discretization itself, 
we take into account additional parameters like the position of the bound- 
ary, lack of resolution of the grid, etc. In [19], for example, it was shown 
that deplacing the boundary helps to correct numerical error resulting in a 
wrong wave velocity and this displacement has been directly introduced in 
the discretization of derivatives. 

In this paper, we apply 4D-Var data assimilation to control the discretiza- 
tion of derivatives and interpolation operators in the boundary regions of a 
non-linear shallow-water model. We use methods proposed in [20] and study 
the assimilation results obtained with a model that exhibits a chaotic be- 
havior. Two examples are considered in this paper: an academic case of 
the model in a square box with artificially generated observations and more 
realistic case of assimilation of real observational data in the Black sea model. 
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2 Shallow Water Model 



2.1 Model's equations and discretization 

The data are assimilated into the shallow-water model on the /^-plane 



dt dy dx PqHq 

drj drju dr]v ^ 
dt dx dy 

where u{x,y,t) and v{x,y,t) are two velocity components, ri{x,y,t) is the 
sea surface elevation, po is the mean density of water, Hq is the characteristic 
depth of the basin and g is the reduced gravity. The model is driven by the 
surface wind stress with components Tx{x,y) and Ty{x,y) and subjected to 
the bottom drag that is parametrized by linear terms au, av. Horizontal 
eddy diffusion is represented by harmonic operators fiAu and fiAv. Coriolis 
parameter is supposed to be linear in y coordinate and is presented as (/o + 
Py). The system is defined in some domain Q with characteristic size L 
requiring that both u and v vanish on the whole boundary of Q. No boundary 
conditions is prescribed for rj. Initial conditions are defined for all variables 
u, V and i]. 

We discretize all variables of this equation on the regular Arakawa's C- 

r 

grid [25] with constant grid step /i = -jy in both x and y directions (see 

figHD 

Mij_i/2(t) = u{ih,jh-h/2,t)foTi = 0,...N,j = 0,...,N + l 

Vt-i/2,j{t) = v{ih - h/2, jh,t) ior i = 0,...N + = 0,...N 

m-i/2,3~-i/2{t) = r]{ih-h/2,jh-h/2,t) hi i = 0,...N + =0,...N + 1 

In order to discretize the system ([T]), we replace the derivatives by their 
finite difference representations and Dy and introduce two interpolations 
in X and y coordinates and Sy. Interpolations are necessary on the stag- 
gered grid to calculate the variable's values in nodes where other variables 
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Figure 1: Arakawa C-grid 
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are defined. Following we write the discretized system ([T]) as 



du T 
-Trr-{h + Py-SyDyu)S,SyV + D,{{Sxuf/2 + gr]) = -au + iiA^u+ " 



^ + ifo + Py + SxDxv)SySxU + Dy{{Syv)y2 + gr]) = -av + fiA^v + 

^ + DxiuSxT]) + DyivSyT]) = 0. 

Discretized operators D^jDy and Sx,Sy are defined in a classical way 
at all internal points of the domain. For example, the derivative and the 
interpolation operator of the variable u defined at corresponding points write 



,Dxu)i-i/2,j-i/2 = — for « = 2, . . . , iV - 1, 



{OxU)i^i/2,j-i/2 — tor t — z, . . . ,I\ - i. 



(3) 



These expressions represent a well known second order approximation. Laplace 

Vi- 

1? 



operator is discretized in a common way A^'i; = — + 



J? • 

The boundary region is considered as a band of adjacent to boundary 
nodes. Discretizations of operators in this band are different from ([3]) and 
represent the control variables in this study. In order to obtain their optimal 
values assimilating external data, we suppose nothing about derivatives near 
the boundary and write them in a general form 

{IJxU)l/2,j-l/2 = (4j 

This formula represents a linear combination of values of u at two points ad- 
jacent to the boundary with coefficients a. These coefficients are considered 
as particular for each operator and for each variable. The constant ao, which 
is also particular for each operator, may be added in some cases to simulate 
non uniform boundary conditions like u{0,y) = a^""" 7^ 0. 

We distinguish a for different variables and different operators allowing 
different controls of derivatives because of the different nature of these vari- 
ables and different boundary conditions prescribed for them. It is obvious. 
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for example, that the approximation of the derivative of rj in the gradient 
may differ from the approximation of the derivative of u in the divergence. 
Although both operators represent a derivative, boundary conditions for u 
and rj are different, these derivatives are defined at different points, at differ- 
ent distance from the boundary. Consequently, it is reasonable to let them be 
controlled separately and to assume that their optimal approximation may 
be different with distinct coefficients and a^""^. 

Time stepping of this model is performed by the leap-frog scheme for all 
hyperbolic terms and Euler scheme for the dissipative terms. The first time 
step is splitted into two Runge-Kutta stages in order to ensure the second 
order approximation. 



2.2 Tangent and adjoint models 

The approximation of the derivative introduced by ([3]) and (jl]) depends on 
control variables a. Operators are allowed to change their properties near 
boundaries in order to find the best fit with requirements of the model and 
data. To assign variables a we shall perform data assimilation procedure and 
find their optimal values. Variational data assimilation is usually performed 
by minimization of the specially introduced cost function. The minimization 
is achieved using the gradient of the cost function that is usually determined 

by the run of the adjoint to the tangent linear model. 

Developing the tangent and adjoint models in this case, we follow the 
procedure presented in [20]. The tangent model is equal to the Gateaux 
derivative of the original model ([2]) with respect to the control parameters. 
To calculate this derivative we suppose that the control variables a can have 
small variations Sa and we determine the linear model that governs the 
evolution of the perturbations 6u, Sv, 6r] of the solution u, v, rj of the 
model ([2]). Skipping the detailed development of the tangent model (see 
[20]). we write the model as 



d_ 

di 



5u 
Sv 
5ri 



where 



SxSyV ■ SyDy 

UO + Py + S^D^V) ■ SyS:, 





— {/O + Py — SyDyU) ■ SxSy 

3y Sx W ■ Sx Dx 
+ Dy{SyV ■ Sy) - a + ^J,A'^ 
Dy{Syr)-) 





-gDx Rr, 
-gDy Rv 








(5) 



Ru5a = -(/O + - SyDyU) ■ {5Sx{SyV) + Sx{5SyV)) + S^^SyV ■ {5Sy{DyU) + Sy{5DyU)) + 

+ 5D,{{S^uf/2 + gTi) + D^{S.,u-5S,u) 
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RJa = {fo + f3y + S^D^v) ■ {SSy{S^u) + Sy{SS^u)) + SyS^u- {SSx{D^v) + S^{SD 

+ 6Dy{{SyVy/2 + gr]) + Dy{Syv6Syv)) 
R^5a = 6D^{u ■ Sa^f]) + Dx{u ■ dSj^rf) + 5Dy{v ■ Syi]) + Dy{v ■ SSyi]) 

Operators 6S and 6D denote S{a + 6a) — S{a) and D{a + 6a) — D{a) 
respectively. These operators are of implicit structure. Their argument (that 
is, in fact, 6a) is contained in the matrix itself. This representation is not 
convenient in the development of the adjoint model, that's why we would 
better rewrite them in a more classical form: a constant operator (which 
does not depend on 6a) multiplied by a variable vector (which is 6a). So, 
each product like 6DxU, 63^^] etc. is replaced by another product: 6DxU = 
[u]^^6a, 63x1] = [ri]g^6a etc. where operators [u], . . . are constructed 

from the solution u, t] of the original equation and the vector 6a is extracted 
from matrices 6D or 63 as it is described in [12]. We shall further use hats 
to denote matrices that have been constructed from vectors. These matrices 
are also block-matrices. All elements of their blocks are equal to zero except 
one line in the beginning and one line at the end of each block. These lines 
are composed of two values of corresponding vector, and, namely, values of 
approximated function in two nodes near the boundary. 
Using these notations, operators R are rewritten as 

Ru = -ifo + /3y - 3yDyu) ■ {[3yv]g^ + 3x[v]s^) + 3x3yV {[Dyu]s^ + 3y[u]^^) + 

+ [{3xu)y2 + gr^]^^ + D.(S.w ■ 

Rv = {fo + Py + 3xDxv)-{[3xu]s^ + 3y[u]gJ + 3y3xU-{[Dxv]g^ + 3x[v]^J + 

+ [{3yv)y2 + g^]^^+Dy{3yv\^],J (7) 

Rrj = [u- 3xri]jy,^+ Dxiu-[ri]gJ + [v 3yri]^^ + Dy{v[ri]gJ 

The tangent model is subjected to the same zero boundary conditions 
for 6u and 6v as for u and v. No boundary conditions is prescribed for 6r) 
as well as for t]. At initial time 6u, 6v and 61] are equal to 6uo, 6vq and 6t]q 
respectively. These variables are now classical in controlling initial conditions 
of the model. They will be used for the same purpose and for the joint control 
of a boundary scheme and initial conditions of the model as well. 

We should note, the matrix of the tangent linear model ([5]) is composed 
by two parts: the 3x3 block composed of operators acting in the space of 
the model's variables and the fourth column composed of operators R ([7]). 
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The 3x3 block is responsible for the evolution of a small perturbation by 
the model's dynamics and is similar for any data assimilation. The column 
determines the way how this perturbation is introduced into the model and 
is specific to the particular variable under control. This column is absent 
when the goal is to identify the initial conditions of the model because the 
uncertainty in initial conditions is introduced only once, at the beginning of 
the model integration. But, when the uncertainty is presented in an internal 
model parameter, like in this case, the perturbation is introduced at each 
time step of the model. 

In order to develop the adjoint model, we need to introduce the scalar 
product in the space defined by tangent model. Each element in this space 
is composed of discretized functions u, v and rj and also the whole set of 
the control coefficients a. A vector in this space has four components (p = 
{u,v,'f],a). 

Following [20], we consider a weighted Euclidian scalar product in this 
space 



/ u \ 




/ u* 


\ 


V 




V* 




7] 




T]* 




V a J 




\ a* 


J 



a 



operator / operator \ * 



hJ 



hJ 



k, operator 



The sums in the scalar product is performed over all nodes i,j of the grid 
of all model's variables u, v and rj. The sum of control coefficients a is 
performed over all : < A; < 2, (|4]) and over all operators ^^operator" 

controlled by these coefficients. Weights Wu = = and w„ = -n- 

are introduced to bring all variables to dimensionless form. Applying the 
development described in jJO] we get the formulation of the adjoint model: 



d_ 

"at 



where 

Rl = 



( 



D;s;{s^Syv)+ 

+S*[S-^u-D%)-a + i,A'^ 
-S;S*{fo + I3y - SyDyu)- 

R*. 



SiS-*{fo + I3y + S^Da:vy 

D*S*{SyS:,U-) + 

+s;{Syv ■ d;) - a + i^Af^ 

H,, 



S.V-D* \ 

Syrj ■ d; 





R* 



1/ 



(9) 



-i[Syv]s^ + [v]s,S:)iifo + Py- SyDyu)-) + {[Dyu],^ + [u]^^S;){S^Syv) + 
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R 



+ [{S,u)y2 + gr^]^^ + [u]s^{S,,u-D:) 



2.3 Cost function 

One of the principal purposes of variational data assimilation consists in the 
variation of control parameters in order to bring the model's solution closer 
to the observational data. This implies the necessity to measure the distance 
between the trajectory of the model and data. Introducing the cost function, 
we define this measure. Generally speaking, the cost function is represented 
by some norm of the difference between model's solutions and observations, 
eventually accompanied by some regularization term. 

To characterize the difference between the model's solution and the ob- 
servational data we use the norm based on the scalar product with the 
fourth component a equal to zero: 



obs\2 



Expressing ^ in terms of the scalar product, we emphasize its dependence on 
time and on control coefficients a: 



e 



^\a, t) =^(j){a, t) - (f)°'''{t), (f){a, t) - </.'"'*(t)»= 

^obs^-l-\ \ / „,fn, +\ r,,obs( 
,obs I 



r/(a,t)-r/«^^(t) 




; 



(12) 



Taking into account the results obtained in [20], we define the cost function 
as 

T 

X{a) = J t^'^{a,t)dt (13) 



that gives bigger importance to the difference at the end of assimilation 
interval. 
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It should be noted here, that this cost function can only be used in the case 
of assimilation of a perfect artificially generated data. When we assimilate 
some kind of real data that contains errors of measurements and is defined on 
a different grid, we should add some regularization term to the cost function 
(like the distance from the initial guess) and use some more appropriate norm 
instead of the Euclidean one (see, for example for details). 

To calculate the gradient of the cost function, we calculate first its vari- 
ation: 

T 

51 = I{a + 6a)-I{a) = Jt{^'^{a + 6a,t)-^'^{a,t))dt = 



T 

= J t(^<^(j){a + Sa,t) - (j)°^'{t),(f){a + 5a,t) - (f)°'''{t)^ - 



T 

- «c0(a,t) - (j)°'"{t),(j){a,t) - (t)"^'{t):^^dt ^ j t <^6(f){t),(j){a,t) - (f)"'" {t):^li1 



As it has been shown in [20j, the scalar product <^6(f)(t), (j){a, t)— 0"''^(t)» 
is equal to «c50(O), ^*(t, O)(0(a, t) - 0°''"(t))» where ^*(t,0) denotes the 
adjoint model (Q integrated from the state 0(a, t) — 0°**(t) back in time from 
t to 0. So, the result of the adjoint model run, being scalarly multiplied by 
50(0) provides the variation of the cost function 

T 

51 = 2 ^S(f){0), J tA*{t,0){(j){a,t) - (f)'''"{t))dt:$> (15) 



As it has been mentioned above, vector 50(0) is composed of 4 components: 
Suq, 6vo, (5?7o and 6a. If we want to control the boundary scheme only, we 
put zero to the first three components of 50(0). Only the fourth component 
of the variation of the cost function (and its gradient) is different from zero 
in this case and only this component is used in the control. On the other 
hand, if our purpose is to control initial state of the model, then vanishing 
6a must be imposed and the first three components of the gradient must be 
used. And, for the joint control of boundary and initial conditions of the 
model we use all four components of the gradient. 
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Thus, the gradient of the cost function writes 

T 

VX = 2 J tA*{t,0){(f){a,t) - (f)°'''{t))dt. (16) 



This gradient is used in the minimization procedure that is implemented 
in order to find the minimum of the cost function: 

if we control the discretization 
of operators near the boundaijv . 
using only the fourth compo- 
nent of the gradient (fT6!) : 
if we control only the initial 
state of the model using thr^iS) 
components of the gradient; 
if we perform the joint control 
of both the initial state and 
the boundary discretization ()f9) 
the model using all four com- 
ponents of the gradient. 

Coefficients a are considered as coefficients achieving an optimal dis- 
cretization of the model's operators in the boundary regions. We use the min- 
imization procedure developed by Jean Charles Gilbert and Claude Lemarechal, 
INRIA [28]. The procedure uses the limited memory quasi-Newton method. 



X(a) = mmX{uo,vo,r]o,a) 



I{uo,vo,T]o) = mm I{uo,vo,T]o,a) 

uo,vo,rjo 



Z{uo,Vo,r]o,a) = min X{uo,Vo,r]o,a) 

uo,vo,rio,a 



3 Model in a square box. 

We start from the data assimilation in frames of the very well studied "aca- 
demic" configuration. Several experiments have been performed with the 
model in a square box of side length L = 2000 km driven by a steady, zonal 
wind forcing with a classical sinusoidal profile 

27r(y - L/2) 
= To cos 

that leads to the formation of a double gyre circulation p9]. The attractor 
of the model and the bifurcation diagram in a similar configuration has been 
described in [3D]. Following their results, we intentionally chose the model's 



12 



parameters in order to ensure chaotic behavior. The maximal wind tension 

on the surface is taken to be rn = 0.5 '^^^f . The coefficient of Eckman 

cm 

dissipation and the lateral friction coefficient are chosen as a = 5 x lO^^s^^ 

2 

and /i = 200^^ respectively. 

As it has been already noted, the Coriolis parameter is a linear function 
in y with /o = 7 x 10~^s~^ and /3 = 2 x 10~^^(ms)~^. The reduced gravity 
and the depth are respectively equal to g = 0.02^, Hq = 1000m. 

All operators in the model are approximated with the second order accu- 
racy both in the interior of the domain and near its boundary. That means 
the expression (jl]), that is used to interpolate functions and to calculate 
their derivatives near boundary, is written with af = —1/h, = 1/h for 
all derivative operators and af = af = 1/2 for all interpolations. That 
gives, for example, the value of the derivative of u at the point i = 1/2 as 

[UxU)l/2,j-l/2 — Jl • 

The solution of the model in this configuration possesses a boundary layer 

near the Western boundary. This is a well known Munk layer [31] that is 

characterized by the local equilibrium between the /3-effect and the lateral 

dissipation. It's width is expressed by the formula ^?=(/i//3)^/^ = 52.3 km 

v3 

for the given parameters. 

As it has been discussed in [32j, the model must resolve this layer with 
at least one grid point (optimally, more than one grid point) in order to 
maintain numerical stability. The work of [33] emphasized the importance 
of having at least two grid points in the Munk layer in order to minimize the 
level of spurious oscillations visible in the velocity fields as well in the field 
of the sea surface elevation. 

The resolution of the model in this section is intentionally chosen to be 
too coarse to resolve the Munk layer. The model's grid is composed of 30 
nodes in each direction, that means the grid-step is equal to 67 km, that is 
more than the Munk layer's width. As it can be seen in figJH there is only 
one grid node in the layer for variables v and t] and no nodes at all for u 
variable. 

Artificial "observational" data are generated by the same model with all 
the same parameters but with 9 times finer resolution (7.6 km grid step). The 
fine resolution model, having 7 nodes in the Munk layer, resolves explicitly 
the layer and must have no spurious oscillations. All nodes of the coarse 
grid belong to the fine grid, consequently, we do not need to interpolate 



13 



V(x) (cm/s) 



Velocity "v(x)" at y=500 km 



30- 
25- 
20- 
15- 
10- 
5- 

-5 

-10 

-15 



hS67km 
liilekm 




500 1000 1500 2000 



X (km) 



Figure 2: Profiles of tlie velocity f (x, y) at y = 500 km obtained on tlie 
coarse grid witfi h = 67km (solid line), and fine grid with h = 7.6 km 
(dashed hne) 



"observational" data to the coarse grid. We just take values in nodes of the 
high resolution grid that correspond to nodes on the coarse grid. A profile 
of the meridional velocity component v is shown in figE] in order to illustrate 
effects of resolved and under-resolved Munk layer. These effects are the most 
visible in the field of this variable. However, similar oscillations are also 
present in the fields of the velocity u and the sea surface height (SSH) t]. 
As it is shown in figUl the point on the boundary [x = 0) does not belong 
to the grid where v variable is defined. However, boundary conditions are 
prescribed at this point imposing vanishing v. 

If we compare profiles of the velocity v{x,y) plotted in figj2] at y = 500 
km (one quarter of the side length), we see that the profile obtained using 
the model with fine resolution is smooth, while the profile on the coarse 
resolution (solid line in figj2]) shows strong spurious oscillations due to the 
unresolved Munk layer. As expected, one point in the layer is not sufficient 
to suppress the numerical mode in the velocity field. We shall assimilate 
data of the fine resolution model in order to bring the solution of the coarse 
resolution's model closer to fine resolution's one. 

The model on the fine grid has been spun up from the rest state during 
3 years and integrated for the subsequent 3 years. From the result of this 
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integration we have extracted values of all three variables at all grid points 
that belong to the coarse grid (as it has been noted, the grids have been 
chosen so, that all grid points of the coarse grid belong to the fine grid). 
This set is used as artificial observations in the following experiments. 

So far the model is nonlinear with intrinsicly unstable solution, there is no 
hope to obtain close solutions in long time model runs because any difference 
(even infinitesimal) between two models grows exponentially on these time 
scales. Consequently, we have to confine our study to the analysis of relatively 
stable properties of the solution. So, we perform two experiments. The first 
one describes a short time evolution of the model's solution simulating the 
forecasting properties of the model. The second one addresses the climatic 
averages of the model's solution which should also represent more stable 
structures than particular trajectories (see [31], [35] for example). 

In the first experiment we analyze the data assimilation that was used to 
identify optimal initial conditions of the model uo,vo,f]Q, optimal discretiza- 
tion of its operators near the boundary a or both applying minimization 
procedure and controlling different parameters as it is shown in f lT7|l .f ll8p 
and 0191] . The final point of spin- up of the high- resolution model was used as 
initial guess in experiments that control the initial state of the low-resolution 
model. Classical second order discretization of all operators near the bound- 
ary ao = 0, af = —1/h, = 1/h, af = af = 1/2 was assumed as initial 
guess in all experiments that control the discretization. 

In the experiment that analyses the forecasting properties of the model, 
we assimilate "observational data" during 5 days and we examine the evolu- 
tion of the difference between the model's solution and observations during 
next 15 days. 

One can see in figJ2]that when we control initial conditions of the model, 
the control moves the initial state (dashed line) far from initial guess. The 
initial difference with observations becomes ^(0) = 0.1. However, at the 
end of the assimilation window the distance from observations is reduced 
to ^(5 days) = 0.07. After the end of assimilation the distance from ob- 
servations increases rapidly. The line approaches the upper solid line that 
corresponds to the solution of the model with no assimilation at all. This 
fact can easily be understood. The model's dynamics has not been modified 
by the assimilations procedure. The model remains the same, possessing the 
same deficiencies. Consequently, it is not surprising that beyond assimilation 
window the solution tends to the attractor of the coarse resolution model. 
The model starts to develop spurious oscillations shown in fig J2] moving away 
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Figure 3: Evolution of the distance ^(t) ( fTTl) during and after assimilation. 



from the attractor of the high resolution model where these oscillations are 
absent. 

Controlling the discretization of model's operators, we modify the model. 
Assimilation, in this case, brings the model dynamics towards the dynamics 
of the high resolution model (used to produce artificial observations). This 
control does not modify the initial point (.^(0) = 0) but the distance ^(5 days) 
at the end of assimilation time (dotted line in fig|3]) is almost equal to the 
distance obtained controlling the initial state of the model. However, beyond 
the window, we see that the modified model's dynamics allows the solution 
to remain closer to observations. The difference ^ on the twentieth day is 
more than two times lower. If we suppose that we assimilate data in the past 
5 days in order to deliver a forecast for 15 days in the future, we see that 
controlling coefficients a provides two times better result. 

If we control both initial state and boundary discretization, we obtain 
dashed line which is similar to the dotted one. Indeed, the assimilation 
procedure was particularly concentrated on the control of coefficients a rather 
than on the control of the initial state, that remains rather close to the initial 
guess for this state (^(0) = 0.01) . In this experiment we have lower distance 
^ at the end of the assimilation window, but almost the same ^ at the end of 
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the forecasting time. 

Thus, we see in this experiment that if the model's dynamics suffers 
from low resolution and other numerical errors, better forecast is achieved 
by controlling the model's operators rather than initial conditions. 

If we consider a long time behavior of the model, we should analyze 
climatic averages of the solution rather than the difference between particular 
trajectories. Due to intrinsic instability any trajectories diverge and the 
value of ^ becomes determined by the attractor size. If we assimilate data 
in order to control initial state of the model, we can not hope to improve 
its climatic averages because we do not modify the model's dynamics. No 
matter from which state the model starts, the same dynamics determines the 
same attractor and the same climatic (calculates over sufficiently long time 
interval) averages. In the same time, data assimilation performed with the 
purpose to control the discretization of operators near the boundary, does 
modify the dynamics. In this case, together with the short time behavior, 
we can hope to improve the model's climate. 

To study the modification of the climatic averages, we perform another 
experiment starting from the same initial guess, controlling also both initial 
state and parameters a, but with assimilation window T = 100 days. Such 
a large window is necessary to collect an observational information about 
a number of physical processes that determine long-time model behavior. 
Optimal initial point Uo,Vo, f/o and optimal d found in the experiment are used 
in the 1000 days model run that means 10 times the assimilation window. 

Averages of the original coarse resolution model suffer a lot from the 
numerical effects that are present due to insufficient resolution. All fields 
contain spurious oscillations (see figl2]) due to unresolved Munk layer. The 
length of the jet-stream in the middle of the square is about 400 km (800-1000 
km in the high resolution model) and the variability of the model's solution 
is very low. Eddy kinetic and eddy potential energies of the high resolution 
model show approximately 20 times bigger amplitudes. These fields are not 
shown in the paper. 

Optimal discretization of operators in the boundary region allows us to 
correct these numerical errors. Average sea surface elevation fj and eddy 
potential energy EPE 

Vix, y) = ^ r Vix, y, t)dt, EPE = i r{r]{x, y, t) - f]{x, y)fdt (20) 
i JO 1 Jo 

obtained in the "observational" run of the high resolution model and in 
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the run of the low resolution model with optimal boundary discretization are 
presented in fig S] and fig 13 One can see in that both averages and variability 
are very similar. That means, optimal discretization modifies the climate of 
the coarse resolution model bringing it closer to the reference model. 




Figure 4: Thousand days average of the Sea Surface Elevation obtained 
with the high resolution model (left) and with the low resolution model and 
optimal discretization (right). 




Figure 5: Thousand days average Eddy Potential Energy obtained with the 
high resolution model (left) and with the low resolution model and optimal 
discretization (right). 
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4 Model of the Black Sea. 



In this section we use the same model, but all the parameters are defined 
to describe the upper layer circulation of the Black sea. Configuration of 
the model and observational data have been kindly provided by Gennady 
Korotaev from the Marine Hydrophysical Institute, National Academy of 
Sciences of Ukraine, Sevastopol, Ukraine. This configuration is described in 

The model grid counts 141 x 88 nodes that corresponds to the grid box of 
dimension 7860 m and 6950 m in x an y directions respectively. 15 minutes 
time step is used for integration of the model. The Coriolis parameter is equal 
to /o = 10~'^s~^ and /3 = 2 x 10~^^ (ms)~^. Horizontal viscosity is taken as 
fi = 200m^s~^. Using a typical density difference, po — Pi of 3.2kg/m^ , and 
unperturbed layer thickness of Hq = 150m, the Rossby radius of deforma- 
tion is estimated at about 22 km. The grid therefore resolves the mesoscale 
processes reasonably well. 

The model has been forced by the ECMWF wind stress data, available as 
daily averages for the years 1988 through 1999. Dynamical sea level recon- 
structed in |37] was used as observational data in this section. These data 
have been collected in ERS-1 and TOPEX/Poseidon missions and prepro- 
cessed by the NASA Ocean Altimeter Pathfinder Project, Goddard Space 
Flight Center. Observational data are available from the 1st May 1992 until 
1999. These data have been linearly interpolated to the model grid. 

So far the sea surface elevation is the only observational variable available 
in this experiment, we put w„ = = in (fTTj) . Consequently, the differ- 
ence between the model's solution and observations is calculated taking into 
account the variable rj only. 

The behavior of the model solution is not chaotic in this configuration. 
Variability in the model is generated directly by the variability of the wind 
stress on the surface. Consequently, we can compare particular trajectories 
of the model on any time interval because their evolution is stable without 
exponential divergence. 

As it has been already noted, absence of observational data for the velocity 
fields brings us to modify the cost function. We have to add the background 
term in the cost function in order to require the velocity field to be suffi- 
ciently smooth. Otherwise, lack of information about velocity components 
in observational data would result in a spuriously irregular fields obtained in 
assimilation. To ensure necessary regularity of u and v we add the distance 
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from the initial guess to the cost function f ll3p . In order to emphasize the 
requirement of smoothness, this distance is measured as an enstrophy of the 
difference between the initial guess and current state: 



■'smooth 



V dx dy 



where m", f° denote velocity components of the initial guess of the mini- 
mization procedure. Of course, this term is only taken into account in the 
identification of the initial state of the model. 

Moreover, using real observational data requires to add at least one an- 
other term to the cost function. One can see in the Figure 2 in [37], spatially 
averaged sea surface elevation of the Black sea exhibits a well distinguished 
seasonal cycle. That means the mass is not constant during a year, it de- 
creases in autumn and increases in spring. Consequently, if we assimilate 
data during a short time (a season or less), we assimilate also the infor- 
mation about the mass flux speciflc for this season. This flux can not be 
corrected later by the model because the discretization of operators near the 
boundary (that controls the mass evolution) is obtained once for all seasons. 
The mass variation of the Black sea reaches 25 centimeters of the sea surface 
elevation. Assimilating data within one season may, consequently, result in 
a persisting increasing or decreasing of the seal level of order of 50 cm per 
year. To avoid this spurious change of the total mass, we must either take the 
assimilation window of at least one year, or prescribe the mass conservation 
to the model's scheme. One year assimilation window is computationally 
expensive and is not justifled by the model's physics. On the other hand, 
prescribed mass conservation removes just the sinusoidal seasonal variation, 
allowing us to keep all other processes and to choose any assimilation window 
we need. 

To correct the mass flux of the model, we add the following term to the 
cost function 

Imass = I - dt (22) 

Similarly to f l2T|) . this term also ensures the regularity of the solution, but it 
is taken into account when the assimilation is performed for identiflcation of 
the boundary parametrization. It can be noted here that other terms may 
be added to the cost function in order to make a numerical scheme energy 
and/or enstrophy conserving, but we do not use them in this paper. 
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The total cost function in this section is composed of three parts: fll3p . 
and ([22]): 

^total = X + ■yilsmooth + J^^mass (23) 

Coefficients 7 are introduced to weight the information that comes from ob- 
servational data (with X) and an a priori knowledge about mass conservation 
and regularity of the solution. 

This modification of the cost function results, of course, in additional 
terms in the gradient: 

VXtotai = VX+ 271 (D;Dy{u - uo) + DID, {v -Vq)] +2-f2j2 U,j (t) - Vij (0)) • 

(24) 

The model is spun up from the beginning of 1988 to May 1992 using the 
wind tension data on the surface. The state corresponding to the 1st of May 
1992 12h GMT is used as the initial guess in the data assimilation procedure 
controlling initial conditions of the model. The assimilation is performed 
following the procedure f[T7j) with the assimilation window T = 1 day and 
the regularization parameter 71 = 0.04. Such a short window was chosen in 
order to get almost instantaneous state of the model to be used in further 
experiment as an initial state. 

This assimilation provides sufficiently smooth velocity fields (see fig|6] 
left) that contain such a specific features of the Black Sea circulation as 
Western and Eastern gyres, Batumi gyre, filament formation along the Cau- 
casus coast and a formation of the Sebastopol eddy (see [36j for discussion of 
these features). Obtained sea surface elevation field (see figj6] right) is close 
to observational data that was used in the assimilation. The value of the 
difference ^ ffTTj) is less than 0.1, that means the actual average difference 
between observed and reconstructed fields is approximately equal to 30 cm. 

This initial state is used in data assimilation experiments that control 
the discretization of the model's operators near boundary. Of course, short 
assimilation window is not appropriate to identify internal model's parame- 
ters. Observational data must contain information, and especially about long 
time model behavior and must be capable to identify optimal discretization 
which is supposed to be constant in time. We have already shown on the ex- 
ample of the total mass evolution that short assimilation windows may lead 
to a wrong model behavior. However, choosing a long assimilation window 
require much computer time. 
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Figure 6: The model initial state corresponding to the 1st of May 1992 
obtained by data assimilation: velocity (left), sea surface elevation (right). 



In this paper we chose T = 50 days window which is longer than synoptic 
time scales. The minimization of the cost function has been accompanied by 
the mass preserving correction f l22|) with 72 = 0.01. 

As well as in the previous section, we perform 3 experiments controlling 
initial conditions, boundary parametrization and both of them. However, 
due to internal stability of the model solution, we can compare particular 
trajectories on long time intervals. Evolution of the difference ^ between the 
model's solution and observations is shown in figJTl As we can see, optimal 
initial point allows the solution to remain close to observations in the as- 
similation window but not beyond the window. As soon as the assimilation 
ends, the solution goes rapidly toward the solution of the original model and 
becomes indistinguishable from it after 200 days. On the other hand, the so- 
lution of the model with optimally discretized operators remains always closer 
to the observational data than the solution of the original model. That means 
modified model's dynamics allows the solution to approach observations. 

Comparing the computational cost of the data assimilation, we must ac- 
knowledge that controlling the boundary is more expensive. First, the adjoint 
model's run requires approximately double computer time comparing with 
the adjoint model for the identification of the initial state. And second, the 
minimization process requires more iterations. In the experiment shown in 
figJTJ 5 iterations are necessary to minimize the cost function controlling the 
initial state of the model, but we need at least 22 iterations to reach the same 
stopping criterion controlling a. 

As it has been discussed in [19], it is difficult to distinguish principal 
modifications that have been made in the numerical scheme by the data 
assimilation. Various kernels are present in the space of a even in a very 
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Figure 7: Evolution of the distance ^(t) ( ITT]) during and after assimilation. 



simple model like one-dimensional wave equation. These kernels make diffi- 
cult the analysis of the assimilation results because optimal coefficients are 
not unique. Obtained set of a represent just one point in the kernel, while the 
same model behavior can be obtained with any other kernel point. Conse- 
quently, in two similar assimilation experiments we can obtain very different 
sets of optimal coefficients but almost the same model behavior. 

In the present non-linear model it is also difficult to analyze coefficients a 
directly. We have performed several experiments assimilating observational 
data during the same seasons with the same assimilation windows but in dif- 
ferent years (1992, 1993, 1995). Assimilation results (not shown here) reveal 
very close values of the cost function obtained in the minimization procedure, 
very similar evolution of the difference ^ in the assimilation window and be- 
yond it (like shown in fig|7]), but very different sets of coefficients a. As well 
as in experiments with simpler models ([I9],[2n]), this non- uniqueness of the 
optimal discretization coefficients has also its roots in the kernel that exists 
in the space of a. Different points in the kernel correspond to different dis- 
cretizations of the model's operators, that result in almost the same model's 
solutions. 

Instead of analyzing the set of obtained coefficients a, we shall see the 
modification of the solution that this set generates, and namely the difference 
between the velocity of the model with classical boundary (the distance of this 
solution from observations is plotted by the solid line in fig|7]) and the velocity 
with optimal boundary discretization (dashed line in fig|7]). This difference 
has been averaged in time over 200 days time interval in order to reveal 
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persistent modifications of tlie flow produced by tlie optimal discretization. 

This average difference of the velocity is presented in figEl We zoom the 
Southern part of the Black sea because it is in this region the difference shows 
the biggest values reaching 15 while in the middle of the sea it rarely 
exceeds 1 




31 31.5 32 32.5 33 33.5 34 34.5 35 35.5 36 36.5 37 37.5 



0.05 

Figure 8: Difference in the velocity field of solutions with classical and opti- 
mal discretizations. 

We can note several principal features of the flow that have been modified 
by boundary conditions. First, we can see a strong current on the boundary. 
The slip condition (vanishing tangential velocity) has been replaced by a per- 
manent current along the boundary. Moreover, impermeability condition has 
also been modified. The flow is now allowed to leave the domain ensuring, 
however, the global mass balance. One can see a strong persistent vortex 
centered at 42.2°N, 32.8°E which southern part crosses the boundary result- 
ing in not only tangential but normal flux also. Similar vortices with lower 
amplitude can also be seen in places where the boundary changes direction. 
Optimal discretization allows the flow to cross the boundary in places where 
the direction change is not smooth. 

Tangential velocity component is amplified in the direct vicinity of the 
boundary. In these nodes we see a strong eastward flow that was forbidden 
by the boundary conditions in the classical formulation of the model. On the 
other hand, the eastward velocity is lower at nodes distanced by several grid 
cells from the boundary. At these nodes we see westward flow in the difference 
of the optimally discretized and classical models. In fact, the control of the 
discretization of operators near the boundary results in the same phenomenon 



24 



as we have seen above in the experiments with the square box showing the 
Munk boundary layer figI21 A strong tangential current in the boundary layer 
is moved on the boundary allowing more optimal representation of a thin 
current on a coarse grid that brings the model solution towards observations. 

5 Conclusion 

This paper is an extension of the study presented in [20]. Both papers discuss 
the variational data assimilation procedure applied for identification of the 
optimal parametrization of derivatives and interpolation operators near the 
boundary, but now we work with nonlinear models assimilating both artificial 
and real observational data. 

Contrary to linear models, even assimilating artificially generated data in 
twin experiments, we can not obtain a solution that is indistinguishable from 
assimilated data. Due to intrinsic instability, even initially close solutions of 
a non-linear chaotic model diverge. Hence, any difference between models 
in twin experiments grows with time and prevent the solutions to remain 
close. Thus, dealing with nonlinear models, we can control the discretization 
near the boundary. This allows us to correct the numerical errors due to 
insufficient resolution (see fig 12]), to bring the model solution closer to ob- 
servations (see figini figlZ]) and to improve statistical averages of a chaotic 
model's solution (see figJH figIS]). 

Controlling the discretization of operators of the model possesses another 
advantage. The solution with optimal discretization remains closer to obser- 
vational data after the assimilation end than the solution with optimal initial 
conditions. This fact can be observed both in the square box and in the Black 
sea. Starting from the optimal initial point, the trajectory remains close to 
observations in the assimilation window, but the distance with observations 
increases rapidly beyond the window. The solution with optimal discretiza- 
tion remains closer to observations even after the end of assimilation time. 
That means, observational data in the past are more efficiently used for fore- 
casting if we assimilate them to control internal model's parameters than to 
control the initial point. 

Acknowledgments. Author thanks Gennady Korotaev from Marine 
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Sea. 
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